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A HIGH ORDER HDG METHOD FOR CURVED-INTERFACE 
PROBLEMS VIA APPROXIMATIONS FROM STRAIGHT 
TRIANGULATIONS 

WEIFENG QIU *, MANUEL SOLANO t, AND PATRICK VEGA ^ 

Abstract. We generalize the technique of [Solving Dirichlet boundary-value problems on eurved 
domains by extensions from subdomains, SIAM J. Sci. Comput. 34, pp. A497-A5I9 (2012)] to 
elliptic problems with mixed boundary conditions and elliptic interface problems involving a non- 
polygonal interface. We study first the treatment of the Neumann boundary data since it is crucial 
to understand the applicability of the technique to curved interfaces. We provide numerical results 
showing that, in order to obtain optimal high order convergence, it is desirable to construct the 
computational domain by interpolating the boundary/interface using piecewise linear segments. In 
this case the distance of the computational domain to the exact boundary is only 0{h?). 
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1. Introduction. In this paper we present a technique to numerically solve 
second order elliptic problems in domains Vt which are not necessarily polygonal. In 
addition, we deal with domains divided in two regions by a curved interface S. In 
particular we use a high order hybridizable discontinuous Galerkin method (HDG) 
[am where the computational domain do not exactly fit the curved boundary or 
interface. The main motivation of this technique is being able to use high order 
polynomial approximations and keep high order accuracy using triangular meshes 
having only straight elements. 

One of the first ideas in this direction was introduced by [5] for the one-dimensional 
case and then extended to higher space dimensions for pure diffusion [S E] and 
convection-diffusion [9] equations. In their work, the mesh does not fit the domain 
and the distance between the computational domain and the boundary T := dVt is of 
only order 0{h), making this method attractive from a computational point of view. 
In addition, [7] applied this method to couple boundary element and HDG methods to 
solve exterior diffusion problems. However, only Dirichlet boundary value problems 
have been considered since Neumann data can not be handled in the same way as we 
will explain below. We will see that for the Neumann boundary case the proposed 
technique works properly if the computational domain is order 0{h^) away from the 
actual boundary. 

The work presented here focuses first on the treatment of part of the boundary 
where a Neumann data is prescribed. It is important to understand this situation in 
order to extend the ideas to problems having a curved interface. In fact, the trans¬ 
mission conditions at the interface involve jumps of the scalar variable and jumps of 
the normal component of the flux. The computational jump of the scalar variable can 
be treated considering the transferring technique of [9] and the computational jump 
of the normal component of the flux can be handled using the extrapolation method 
for the Neumann data that we will describe in the following sections. 
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One of the first methods that approximate Neumann boundary conditions on 
curved domains considering non-fitted meshes was introduced by [1]. Here, a piecewise 
linear finite element method was considered and optimal convergence in the H^-noim 
was shown. In addition, the same authors solved a semi-definite Neumann problem 
on curved domains using a similar technique ([2]). They showed optimal behavior of 
the errors in and L^-norms using again piecewise linear elements. On the other 
hand, higher order approximation finite element methods require to properly fit the 
boundary in order to keep high order accuracy. For instance, isoparametric element 
can be considered (12],[T2]). In the case of elliptic interface problems, usually the curve 
describing the interface is interpolated by a piecewise linear computational interface. 
Hence, super-parametric elements near the interface must be considered in order to 
achieve high order accuracy m)- 

This article aims to develop a high order method based on a triangulation of the 
domain involving only straight elements. As we will discuss, the boundary/interface 
must be interpolated by piecewise linear function in order to obtain the expected 
rates of convergence. Since most of the methods based on linear fitting are only 
second order accurate, we believe our method constitutes a competitive alternative. 

The rest of the manuscript is organized as follows. We will begin by setting 
notation. Then, we will describe the technique for a boundary-value problem where 
Neumann data is prescribed in part of the boundary. In particular, we will discuss the 
proper choice of the paths that will transfer the Dirichlet and impose the Neumann 
data. We will provide numerical simulations showing the performance of the method. 
Then, we will adapt these ideas in order to solve a elliptic interface problem and show 
numerical experiments validating the technique. 


2. Mesh construction and notation. Let be a a triangulation constructed 
by the union of disjoint straight triangles that approximates a bounded domain Q C 
M? and does not necessarily fit its boundary. The Dirichlet and Neumann part of the 
boundary T are denoted by T^:) and (ri:)nr 7 v = 0, ri:)Ur 7 v = r). We also assume 
that the computational boundary, T^, satisfies T^ = T^ UT^ and T^ flT^ = 0 where 
and are part of T^ with Dirichlet {qd) and Neumman (qn) data, respectively. 
Let (i(r,r^) be the distance between T and T^. We denote by Hk the diameter of 
the element K G and by n its outward unit normal. The meshsize h is defined 
as maxxGD^ Let be the set of interior edges of and the edges at the 
boundary. We say that an edge e G if there are two elements and K~ in Dh 
such that e = dK^ D dK~. Also, we say that e G if there is an element K G 
such that e = dKOT^. We set 8/^ = For each element K in the triangulation 

D/j,, we denote by the space of polynomials of degree at most k defined on the 

element K. For each edge e in Eh T^(e) is the space of polynomials of degree at 
most k defined on the edge e. Given an element K, ('r)K and (*, •)aK denote the 
LhK) = {v:f^v^< 00 } and L^(dK) = U : fax < oo} products, respectively. 
Thus, for each f and we define 


^ and (C,'0)aD^ = ^ 

KeDh KeDh 
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3. Boundary value problem with mixed boundary conditions. We con¬ 
sider the following model problem: 


(3.1a) 

(3.1b) 

(3.1c) 

(3.1d) 


- V ■ q = f in n, 
q + K\/u = 0 in 

u = Qd on F/:), 
q n = qn on Fat. 


Here qd ^ and qn ^ H ^/^(Fat) are given data at the border, / G 

is a source term and K G [L ^is a symmetric and positive definite tensor. 


In the computational domain D^, problem (3.1) can be written as follows: 


(3.2a) 

- V • g = / in D^, 

(3.2b) 

g + KVm = 0 in D^, 

(3.2c) 

u = gD on r' 

(3.2d) 

q ■ n = gN on F' 




Here gn and are unknowns. As we mentioned before, gn can be calculated 
following [5l|7l[9], i-c., 


(3.3) 


goix) := goix) + / 

J <j{x) 


K ^q mds^ 


where cr(cc), is a path starting at cc G F^ and ending at x G F^^); and m is the tangent 
vector to (j(x). This expression comes from integrating (3.1b) along the path (j(x) 
(see [9] for details). 

In principle, any kind of numerical method using polygonal domains can be used 
to solve the equations in D^. However, it is desirable to consider those methods 


where an accurate approximation of q is obtained, since the boundary condition (3.3) 


depends on that flux. We also notice from (3.3) that the same idea will not work for g^ 


since a similar expression will involve derivatives of q which are not well approximated 
by the numerical method. 

3.1. The HDG method. The method seeks an approximation (q^^Uh^Uh) of 
the exact solution {q,u^u\^^) in the space x Wh x given by 


(3.4a) 

(3.4b) 

(3.4c) 


Vh = {ve [L\Dh)]^ 
Wh = {we L\Dh) : 
Mh = {fie L^{Eh) : 


v\k e [y\K)f 
w\Key'^{K) 
mU e 3’''(e) 


Vi^ G D J, 

Vi^ G D4, 
Ve G Eh}- 


It is defined by requiring that it satisfies the equations 


(3.5a) - (K-q^, Vw)d^ + {q^ ■ n, w)aD^ = (/, w)d,, 

(3.5b) (g^j, v)o^ - {uh, V ■ v)o^ + {uu, v ■ n)dD^ = 0, 

(3.5c) = 0 , 

(3.5d) {fi,Uh)Yii^ ~ 

(3-5e) (/i, g^)ph = (M!fliv)r'f,; 
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for all {v^w^ii) G Vh x Wh x Here is the approximation of go proposed by 
[9]. More precisely, let K We define the operator ^ 

such that E^{v) = v for all v G Then, for cc G e C T^, 


(3.5f) 


goix) ^ g^{x) := goix) E K E -{q^)-mds, 

J cr(x) 


where is the triangle where e belongs. In other words, E^^ is the standard 
extension of a polynomial to the whole space. On the other hand, g^ is an 
approximation of g^ which is still unknown. In Subsection 3.3 we propose to replace 


(3.5e) by an equation involving known quantities at the right hand side. 


Finally, to complete the definition of the HDG method we must specify the defi¬ 
nition of numerical trace on which we takes of the form 


(3.5g) 


Qh = Qh + Ti^h - Uh)n, 


where r : dD^ (0, oo) is a stabilization parameter that guaranties solvability of 
(3.5) and can be set as t\k = \\^\\l°^{k) oii each element K ( j3l [T3] ). 


3.2. Definition of the family of paths. The representation of g^ in (3.5f) is 
independent on the integration path. Let cc be a point on a boundary edge e. Previous 
work have proposed two ways to determine a point x in T and hence construct (j{x)\ 
(PI) If X is a vertex, an algorithm developed by [9] uniquely determines x as 
the closest point to x such that cr(cc) does not intersect another path before 
terminating at T and does not intersect the interior of the domain O. In 
addition, if x is not a vertex, its corresponding path is defined as convex 
combination of those paths associated to the vertices of e. For the Dirichlet 
boundary value problem, the authors in [9] numerically showed optimal rates 
of convergence with this choice of (7{x) when (i(r,r^) is of order h, that is, 
order /c + 1 for Uh and and order /c + 2 for the numerical trace Uh • 

(P2) On the other hand, [8] proposed to determine x such that m is normal to the 
edge e. In this case these authors theoretically proved that if d(r,r^) is of 
order h, the order of convergence for and is indeed /c +1, but the order 
for Uh is only /c + 3/2. However, if d(r, F^) is of order the numerical trace 
also super converges with order /c + 2. Moreover, they also showed numerical 
evidence indicating that the numerical trace optimally super converges even 
though d(r,r^) is of order h. 


Let now be e a boundary edge with vertices Xi and X 2 - We denote by Fg the 
part of F determined hy Xi and X2 as it is shown in Fig. In this paper we assume 
that if e C F^ (or e C F^) then Te cTd (Fg C Fat). The algorithm in (PI) can be 
easily modified to satisfy this assumption. On the other hand, the paths defined in 
(P2) will not always satisfy this condition. 

3.3. Approximation of the Neumann boundary condition. Let e C F^ a 

Neumann boundary face and Fg C Fat the part of Fat associated to e. We denote by 
Ke the element of the triangulation where e belongs. 

The main idea is to characterize Fg using the parameterization induced by the 
family of paths. More precisely. Let e = {x \ x{0) = {x 2 — Xi)6> + Xi, G [0,1]}. Then 

Te = {x = 4>{e): 4>{e) = x{e) + \a{x{e))\m{e),0 e [o, i]}, 


(3.6) 
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Fig. 1. Examples of a boundary edge e with vertiees xi and X2. Fe is the segment of Vn 
determined by xi and X2- 


where we recall that \a{x{0))\ and m{0) are the length and tangent vector of the 
segment joining x{0) and x{0). We define the space 


(3.7) M^(re) := I M e i^^(re) : M = withG ^^([0,1]) 1. 

Il<^ II2 


Equation (3.5e) is then replaced by imposing the following condition over q^: 
( 3 . 8 ) {E^‘‘{q^) ■ n,fi)r^ = {gN,f^)re V/x G M<^(re). 

Notice that (|3.8|) becomes 


(3.9) [ {q^) ■ n) o 0) {6) fl{0)d0 = [ {qn o 4>) g{0)d0 

JO JO 

for all ft{{0)) G ^/c([0,1]); hence, there is no need of computing the derivative of 0. 


On the other hand, we observe that if m and a were independent of 0 (for example, 
if Eg were polygonal and m perpendicular to e), then ||(/)' o 0“^ ||2 would be constant 
and hence M(f){Te) becomes a standard space of polynomials through pulling back 
polynomials from the interval [0,1]. As we will see in the numerical experiments 
provided in next section, this technique performs optimally if m and n have the same 
direction. 


3.4. Numerical results: boundary-value problem. In this section we present 
numerical experiments showing the performance the extrapolation technique and the 
influence of the choice of paths. Since the size of the computational domain changes 
with h, we measure the errors := u — Uh^ Cq q — and eq := u — Uh using 
the following norms: 


||e«||L2(Dfc) 


ID, 


, 1/2 




l|eg||[L^(Dfe)]=^ 

|D,.|V2 




1/2 
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Here Pa is the projection over IP^(e) with e C dK. 

In addition we compute an element-by-element postprocessing, denoted by of 
the approximate solution Uh, which provides a better approximation for the scalar 
variable when k > i(iaE]). Given an element K we construct ^ as the 

only function in such that 

Uh= I 

I '^hdx if > 0, 

and Uh is the polynomial in (set of functions in with mean zero) 

satisfying 


{yuh,vw)K = -{Qh,^w)K ywey'^+\K). 

In the purely diffusive case, this new approximation of u has been proven to converge 
with order k 2 ioi k > 1 when the domain is polygonal mm), and also when it 
has curved Dirichlet boundary mm)- 


We set K = I in all the experiments of this section. In Subsection |3.4.1| we show 
that deteri orat e convergence can happen if d(r,r^) = 0{h). However, we will see in 
Subsection 3.5 that optimal convergence is obtained when (i(r,r^) = 0{K^). 

3.4.1. Computational domain at a distance (i(r,r^) = 0{h), In the fol¬ 
lowing examples the computational domain is constructed in such a way that the 
distance d(r,r^) is of order h. Moreover, /, go and gjsf are chosen in order that 
u{x,y) = sin(x) sin(^) is solution the exact of (3.1). 


Example 3.1. Our first example consist of approximating a squared domain 
ft = (0,1) by a squared subdomain satisfying (i(r,r^) = 0{h) as Fig. shows. Let 
Tjsf = {x : X = 0}, Td = dfl \Tjsf and the family of paths is computed according to 

(P2). 



_ Fig. 2. Two consecutive meshes (h = l/4: and h = l/S) approximating the domain of Example 

EH (Figure obtained from m 

In Table we display the history of convergence for different polynomial degree 
{k = 0,1,2 and 3) and meshsizes {h = 1/2,1/4,1/8,1/16 and 1/32). We observe 
that the error of u and q behaves optimally with convergence rate of order k 1. 
Moreover the error of numerical trace and postprocessed solution also converge with 
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order /c + 1, which is not optimal for the standard HDG method on polygonal do¬ 
mains. Even though, the errors Cu* are always small than We attribute this lack 
of superconvergence to the fact that the Neumann condition (3.8) is being imposed 
on and not on as in the standard HDG method. 




lle.llint 

Ikq II int 

IkdU 

h 

llenHIi 

nt 

k 

h 

error 

order 

error 

order 

error 

order 

error 

order 


1/2 

4.58E-03 

- 

6.59E-02 

- 

2.13E-02 

- 

7.50E-03 

- 


1/4 

6.09E-03 

-0.41 

4.77E-02 

0.46 

5.75E-03 

1.89 

6.60E-03 

0.18 

0 

1/8 

4.62E-03 

0.40 

2.74E-02 

0.80 

1.75E-03 

1.71 

4.71E-03 

0.49 


1/16 

2.78E-03 

0.73 

1.46E-02 

0.91 

6.18E-04 

1.51 

2.80E-03 

0.75 


1/32 

1.52E-03 

0.87 

7.52E-03 

0.96 

2.51E-04 

1.30 

1.53E-03 

0.88 


1/2 

1.54E-03 

- 

9.89E-03 

- 

3.70E-03 

- 

1.67E-03 

- 


1/4 

5.67E-04 

1.44 

2.55E-03 

1.96 

6.31E-04 

2.55 

4.68E-04 

1.84 

1 

1/8 

1.69E-04 

1.75 

7.09E-04 

1.85 

1.50E-04 

2.07 

1.31E-04 

1.83 


1/16 

4.62E-05 

1.86 

1.94E-04 

1.87 

3.84E-05 

1.97 

3.60E-05 

1.87 


1/32 

1.21E-05 

1.93 

5.13E-05 

1.92 

9.83E-06 

1.97 

9.52E-06 

1.92 


1/2 

2.29E-04 

- 

1.20E-03 

- 

5.23E-04 

- 

2.17E-04 

- 


1/4 

2.82E-05 

3.02 

1.24E-04 

3.28 

3.36E-05 

3.96 

2.44E-05 

3.16 

2 

1/8 

3.43E-06 

3.03 

1.36E-05 

3.19 

3.22E-06 

3.38 

2.81E-06 

3.12 


1/16 

4.25E-07 

3.01 

1.63E-06 

3.06 

3.61E-07 

3.16 

3.38E-07 

3.05 


1/32 

5.28E-08 

3.01 

2.02E-07 

3.01 

4.26E-08 

3.08 

4.13E-08 

3.03 


1/2 

3.37E-05 

- 

1.51E-04 

- 

7.55E-05 

- 

3.39E-05 

- 


1/4 

2.30E-06 

3.87 

9.32E-06 

4.02 

3.12E-06 

4.59 

2.30E-06 

3.88 

3 

1/8 

1.55E-07 

3.89 

6.74E-07 

3.79 

1.78E-07 

4.14 

1.55E-07 

3.89 


1/16 

1.05E-08 

3.89 

4.76E-08 

3.82 

1.12E-08 

3.99 

1.05E-08 

3.89 


1/32 

6.90E-10 

3.92 

3.22E-09 

3.89 

7.13E-10 

3.97 

6.90E-010 

3.92 


Table 1 

History of convergence of the approximation in Example \3.1\ 


Example 3.2. We now consider an annular domain Cl = {{x^y) G : 14^ < 
< 20^} that is being approximated by a polygonal subdomain satisfying 
(i(r,r^) = 0{h) as shown in Eig. i We consider Neumman data in the outer 
boundary Etv = G : x‘^ f- y‘^ = 20^} and Dirichlet data in the inner 

boundary Vd = {{x^y) : = 14^}. Here the paths are computed according to 

(P2). 



Fig. 3. Annular domain and mesh in Example \3. 2\ 
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The behavior of the L^-norm of the error displayed in Table is similar to the one 
obtained in the previous example, i.e., the rate of convergence of the error in all the 
variables is of order k 1. Thus, this example suggests that our technique performs 
properly when the boundary is actually non-polygonal. 


||6^6||int IlCqllint HCulle^ || || int 


k 

h 

error 

order 

error 

order 

error 

order 

error 

order 


1.89 

9.56E+00 

- 

8.79E+00 

- 

4.66E-01 

- 

9.80E+00 

- 


0.96 

8.47E+00 

0.18 

5.82E+00 

0.61 

3.72E-01 

0.33 

8.50E+00 

0.21 

0 

0.49 

5.72E+00 

0.57 

3.38E+00 

0.79 

2.42E-01 

0.63 

5.72E+00 

0.56 


0.24 

3.29E+00 

0.81 

1.82E+00 

0.90 

1.37E-01 

0.83 

3.29E+00 

0.81 


0.12 

1.76E+00 

0.91 

9.42E-01 

0.91 

7.26E-02 

0.92 

1.76E+00 

0.91 


1.89 

2.03E+01 

- 

7.85E+00 

- 

9.56E-01 

- 

2.04E+01 

- 


0.96 

5.94E+00 

1.82 

2.12E+00 

1.94 

2.58E-01 

1.94 

5.96E+00 

1.82 

1 

0.49 

1.43E+00 

2.08 

5.03E-01 

2.10 

6.00E-02 

2.13 

1.43E+00 

2.08 


0.24 

3.40E-01 

2.09 

1.20E-01 

2.08 

1.40E-02 

2.11 

3.40E-01 

2.09 


0.12 

8.19E-02 

2.06 

2.92E-02 

2.06 

3.35E-03 

2.11 

8.20E-02 

2.06 


1.89 

4.04E+00 

- 

1.82E+00 

- 

1.90E-01 

- 

4.04E+00 

- 


0.96 

6.80E-01 

2.64 

3.42E-01 

2.46 

2.95E-02 

2.76 

6.81E-01 

2.64 

2 

0.49 

1.41E-01 

2.30 

5.86E-02 

2.58 

5.89E-03 

2.36 

1.41E-01 

2.30 


0.24 

2.12E-02 

2.75 

8.33E-03 

2.83 

8.75E-04 

2.77 

2.12E-02 

2.75 


0.12 

2.88E-03 

2.89 

l.lOE-03 

2.93 

1.16E-04 

2.90 

2.88E-03 

2.93 


1.89 

4.12E+00 

- 

1.52E+00 

- 

1.93E-01 

- 

4.12E+00 

- 


0.96 

3.17E-01 

3.80 

1.07E-01 

3.93 

1.37E-03 

3.92 

3.17E-01 

3.80 

3 

0.49 

1.89E-02 

4.13 

6.29E-03 

4.15 

7.89E-04 

4.18 

1.89E-02 

4.13 


0.24 

l.lOE-03 

4.13 

3.70E-04 

4.12 

4.53E-05 

4.15 

l.lOE-03 

4.13 


0.12 

6.56E-05 

4.08 

2.23E-05 

4.07 

2.68E-06 

4.09 

6.56E-05 

4.08 


Table 2 

History of convergence of the approximation in Example \3.2\ 


Remark 3.1. The construction of the family of paths according to (PI) in Ex¬ 
amples 3.1 and 3.2 deliver similar results since the difference between (PI) and (P2) 
is not significant for these domains. That is why we do not display the convergence 
tables for this case. This numerical evidence indicates that the technique proposed 
provides optimal rate of convergence when d(r,r^) = 0(h) and the family of paths is 
constructed according to (PI) or (P2). However, in practice, this condition over the 
distance can not be satisfied in general, unless the mesh is constructed properly to do 


A practical construction of the computational domain was described in m- It 
consists of ^Immersing the domain in a Cartesian background mesh and set as 
the union of all the elements that are completely inside of ft as it is shown in Fig. 

Here d(T,T^) = 0(h). In this case it is not convenient to construct the paths 
according to (P2). In fact, given a point x e it might happen that x is extremely 
far from x, specially in parts ofT where the domain is non-convex. Since both proce¬ 
dures deliver similar results in previous examples, we will consider from now on (PI). 


Example 3.3. In order to observe the performance of the method where the 
mesh satisfies d(r,r^) = 0(h) and the paths are given by (PI), we consider the ring 
Q = {(x,y) G : 0.25^ < (x — 0.5)^ + (^ — 0.5)^ < 1} with Tn = {(x,y) G : 
x‘^ y‘^ = 1} and To = {(x^y) '. x‘^ -\- y^ = 0.25^}. In Fig. [^we show a zoom at 

the upper-right corner of three consecutive meshes. We also plot the family paths 
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Fig. 4. Left: Domain Cl, its boundary F (solid line), a haekground mesh and the polygonal 
subdomain (gray). Right: Diriehlet data g on V transferred to cp on F/^,. (Figure taken from [P|/) 


from vertices and quadrature points on the boundary edges. In Table we display 
the history of convergence. Even though the method is still convergent for /c = 0, 1 
and 2, the rates deteriorate. Moreover, there is no convergence when k = 3 . For 
the Diriehlet boundary value problem this non-optimal behavior does not occur as 
[9] showed. This example suggests that in a practical situation (meshes satisfying 
d(r,r^) = 0{h) and paths constructed using (PI), the method does not perform 
properly. So, it seems that for Neumann boundary data, the family of paths needs to 
be build according to (P2). Even though we have no theoretical support that explains 
this behavior, we believe it might be related to the oscillatory nature of high degree 
polynomials. In fact, for the Diriehlet problem, [8] showed error estimates where 
some of the constants depend on the polynomial degree. In addition, m numerically 
studied the robustness of this method applied to a convection-diffusion problem with 
Diriehlet boundary data. The concluded that, even though d(r,r^) = 0(h), T and 
must be “close enough’ when k > 1. 

One way of always being able to construct the paths using (P2) is to interpolate 
the boundary by a piecewise linear function. In this case d(r,r^) = 0{h?). 

Remark 3.2. In Example \3.3\ it is not possible to construct the family of path by 
(PI). In fact, a path perpendicular to an inner boundary edge might not intersect the 
inner ring . Moreover, a path perpendicular to an outer boundary edge might intersect 
the outer boundary extremely “/ar’’ as would happen in the third mesh of Fig. 

3.5. Computational domain at a distance d(r,r^) = Another prac¬ 

tical construction of Dh is defining first T^ by interpolating T using piecewise linear 
segments. Then, is the domain enclosed by as Fig. shows. In this case 
(i(r,r/i) = 0(Ji?) and the family of paths can be easily defined according to (P2). 

Example 3.4. We consider the domain Q = {{x,y) G : 1 < (x — 0.5)^ + (^ — 
0.5)^ < 4} with r TV = {{x, y) G : x‘^ y‘^ = 1} and Td = {{x, y) : x‘^ y‘^ = 4}. 

In Table we observe again that the order of convergence in all the variables in 
/c + 1. We point out that part of the computational domain is outside of D as it can 
be observed in the inner circle in Fig. This was never the case in the examples 
provided by [9] and [8]. Thus, these results indicates that their technique also works 
when n D/j, 7 ^ 0. In Fig. we show the approximated solution ph considering 
h = 1.10 (left) and 0.55 (right) and using polynomials of degree /c = 0,1 and 2. We 
clearly see an improvement either when the mesh is refined or the polynomial degree 
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Fig. 5. Zoom at the upper-right eorner of three eonsecutive meshes of Example \3.3\ Mesh (grey 
region) eonstrueted eonsidering the proeedure in (Pf and family of paths determined aeeording to 
(PI). Blue lines: paths from the vertiees. Red lines: paths from quadrature points of the boundary 
edges (k = 1). 


k 

h 

llenlli: 

error 

nt 

order 

\\eq ||i] 

error 

nt 

order 

lleslU 

error 

h 

order 

lle.Hh 

error 

int 

order 


0.312 

4.12E-02 

- 

1.83E-01 

- 

4.40E-02 

- 

4.15E-02 

- 


0.156 

3.70E-02 

0.16 

1.27E-01 

0.53 

3.26E-02 

0.43 

3.69E-02 

0.17 

0 

0.078 

1.69E-02 

1.13 

1.37E-01 

-0.11 

1.50E-02 

1.12 

1.69E-02 

1.13 


0.039 

9.11E-03 

0.89 

7.00E-02 

0.96 

7.61E-03 

0.97 

9.11E-03 

0.89 


0.019 

8.50E-03 

0.10 

4.92E-02 

0.51 

5.66E-03 

0.43 

8.50E-03 

0.10 


0.312 

6.13E-03 

- 

1.82E-02 

- 

3.75E-03 

- 

5.71E-03 

- 


0.156 

3.44E-03 

0.84 

1.06E-02 

0.77 

2.18E-03 

0.78 

3.37E-03 

0.76 

1 

0.078 

3.86E-03 

-0.17 

9.41E-03 

0.18 

2.36E-03 

-0.11 

3.86E-03 

-0.20 


0.039 

1.16E-03 

1.74 

2.68E-03 

1.81 

6.88E-04 

1.78 

1.16E-03 

1.73 


0.019 

5.17E-04 

1.16 

1.16E-03 

1.20 

3.04E-04 

1.18 

5.16E-04 

1.16 


0.312 

4.68E-04 

- 

1.25E-03 

- 

3.03E-04 

- 

4.60E-04 

- 


0.156 

2.25E-04 

1.06 

5.89E-04 

1.08 

1.45E-04 

1.06 

2.24E-04 

1.04 

2 

0.078 

1.21E-04 

0.89 

3.24E-04 

0.86 

7.39E-05 

0.97 

1.21E-04 

0.89 


0.039 

1.31E-05 

3.20 

3.60E-05 

3.17 

7.79E-06 

3.25 

1.31E-05 

3.21 


0.019 

2.63E-06 

2.32 

7.03E-06 

2.35 

1.54E-06 

2.33 

2.63E-06 

2.32 


0.312 

3.02E-05 

- 

8.78E-05 

- 

1.98E-05 

- 

3.00E-05 

- 


0.156 

l.llE-05 

1.44 

3.45E-05 

1.35 

7.19E-06 

1.45 

l.lOE-05 

1.44 

3 

0.078 

1.65E-06 

2.75 

5.37E-06 

2.67 

l.OlE-06 

2.83 

1.65E-06 

2.75 


0.039 

6.69E-06 

- 

1.53E-05 

- 

3.98E-06 

- 

6.70E-06 

- 


0.019 

8.03E-03 

- 

2.26E-02 

- 

4.73E-03 

- 

8.04E-03 

- 


Table 3 _ 

History of eonvergenee of the approximation in Example \3. 3\ 


increases. 
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Fig. 6. Zoom at the upper-right eorner of Example \3.4\ Blue line: boundary F. Grey region: 
mesh. 



l|e„||i: 

nt 

II Cq II int 

llealU 

■h 

IknHh 

int 

k 

h 

error 

order 

error 

order 

error 

order 

error 

order 


1.72 

5.31E-01 

- 

2.14E+00 

- 

2.22E-01 

- 

6.63E-01 

- 


1.10 

2.87E-01 

1.37 

1.19E+00 

1.3 

1.14E-01 

1.48 

3.00E-01 

1.77 

0 

0.55 

1.45E-01 

0.99 

6.13E-01 

0.95 

5.76E-02 

1.00 

1.46E-01 

1.04 


0.29 

8.10E-02 

0.89 

3.31E-01 

0.95 

3.10E-02 

0.95 

8.05E-02 

0.91 


0.15 

4.36E-02 

0.98 

1.69E-01 

1.07 

1.60E-02 

1.05 

4.34E-02 

0.98 


0.08 

2.24E-02 

0.99 

8.48E-02 

1.02 

8.12E-03 

1.01 

2.23E-02 

1.00 


1.72 

2.59E-01 

- 

9.51E-03 

- 

9.51E-03 

- 

1.22E-01 

- 


1.10 

7.11E-02 

2.89 

1.61E-03 

3.97 

1.61E-03 

3.97 

1.80E-02 

4.27 

1 

0.55 

1.77E-02 

2.01 

2.50E-04 

2.68 

2.50E-04 

2.68 

2.54E-03 

2.82 


0.29 

4.45E-03 

2.12 

5.92E-05 

2.22 

5.92E-05 

2.22 

4.23E-04 

2.76 


0.15 

1.08E-03 

2.26 

1.43E-05 

2.25 

1.43E-05 

2.25 

9.03E-05 

2.45 


0.08 

2.66E-04 

2.08 

4.24E-06 

1.81 

4.24E-06 

1.81 

2.69E-05 

1.80 


1.72 

4.59E-02 

- 

6.22E-02 

- 

1.43E-03 

- 

1.04E-02 

- 


1.10 

6.55E-03 

4.35 

9.09E-03 

4.29 

1.95E-04 

4.44 

1.35E-03 

4.56 

2 

0.55 

8.37E-04 

2.97 

1.26E-03 

2.85 

l.lOE-05 

4.15 

8.25E-05 

4.03 


0.29 

1.12E-04 

3.09 

1.71E-04 

3.07 

2.14E-06 

2.52 

1.44E-05 

2.67 


0.15 

1.42E-05 

3.29 

2.11E-05 

3.32 

2.01E-07 

3.75 

1.34E-06 

3.77 


0.08 

1.77E-06 

3.10 

2.63E-06 

3.10 

3.37E-08 

2.66 

2.22E-07 

2.68 


1.72 

5.61E-03 

- 

8.48E-03 

- 

.57E-04 

- 

1.28E-03 

- 


1.10 

4.47E-04 

5.65 

6.59E-04 

5.71 

6.52E-06 

7.11 

4.82E-05 

7.32 

3 

0.55 

3.31E-05 

3.75 

4.77E-05 

3.78 

1.77E-07 

5.20 

1.42E-06 

5.08 


0.29 

2.26E-06 

4.12 

3.30E-06 

4.11 

1.51E-08 

3.78 

1.04E-07 

4.01 


0.15 

1.37E-07 

4.46 

2.12E-07 

4.36 

9.59E-10 

4.39 

6.42E-09 

4.43 


0.08 

8.47E-09 

4.14 

1.32E-08 

4.13 

9.52E-11 

3.43 

6.28E-10 

3.46 


Table 4 

History of eonvergenee of the approximation in Example \3.4\ 


Example 3.5. Now we test the performance of the method where is a bounded 
domain exterior to an airfoil. This is the most difficult case in our examples since the 
domain has a boundary with a curved, re-entrant corner. The airfoil is obtained by 
using the Joukowsky transformation: 

\2 

J{z) = z-\ -, 

z 

where z G C and A G M. It is well known that this transformation maps the disc 
centered at (si,52) of radius R to an airfoil when we set A = i? — \/s\^ Here, 
we take R = 0.1605 and = 52 = 0.01. In Fig. |^we show two triangulations of 
the domain with meshsizes h = 0.143 and 0.073. Neumann boundary conditions are 
imposed around the airfoil and Dirichlet data in the remaining part of the boundary. 
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Fig. 7. Approximation of the scalar variable in Example \3.4\ Columns: meshsize h = 1.10 and 
0.55. Rows: Polynomial of degree = 0,1 and 2. 


We consider the following two examples: 

a) Smooth solution. We set / and g such that u{x, y) = sin(x) sin(^) is the 
exact solution as in previous example. In Table we observe that similar 
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conclusions to those in previous examples can be drawn, even though in the 
case the domain is more complicated. 

b) Non-smooth solution. We now consider a potential flow around the airfoil 

where the exact solution in polar coordinates is u{r,0) = rcos{0) 

Here = 0 around the airfoil. In this case \/u has singularities at the 
leading and trailing edges, hence we do not expect high order convergence 
rates. In fact, this can be seen on Table where in all the cases u converges 
with order one and q converges with order less than one. However, for a fixed 
mesh, the errors decrease when the polynomial degree increases. In Fig. [^we 
show the approximation of the x-component of q considering h = 0.143 and 
0.024 and /c = 0,1 and 2. 




Fig. 8. Meshes of Example \3.5\ Meshsizes h = 0.143 and 0.073. 


k 

h 

error order 

II Cq II int 
error order 

lleslU 

error 

h 

order 

error 

order 


0.143 

5.69E-03 

- 

2.25E-02 

- 

1.35E-03 

- 

5.76E-03 

- 


0.113 

4.78E-03 

0.75 

1.71E-02 

1.18 

7.52E-04 

2.50 

4.81E-03 

0.77 

0 

0.073 

3.12E-03 

0.98 

1.05E-02 

1.11 

4.30E-04 

1.29 

3.14E-03 

0.98 


0.038 

1.59E-03 

1.03 

5.36E-03 

1.04 

1.97E-04 

1.19 

1.59E-03 

1.04 


0.024 

9.93E-04 

1.02 

3.25E-03 

1.08 

1.21E-04 

1.06 

9.94E-04 

1.02 


0.143 

1.41E-04 

- 

2.91E-04 

- 

1.46E-05 

- 

1.48E-05 

- 


0.113 

8.04E-05 

2.38 

1.68E-04 

2.33 

8.36E-06 

2.39 

8.46E-06 

2.37 

1 

0.073 

3.36E-05 

2.01 

6.72E-05 

2.11 

1.95E-06 

3.35 

1.96E-06 

3.36 


0.038 

8.51E-06 

2.11 

1.74E-05 

2.07 

5.30E-07 

2.00 

5.14E-07 

2.05 


0.024 

3.21E-06 

2.11 

6.50E-06 

2.12 

1.32E-07 

3.00 

1.28E-07 

3.00 


0.143 

1.89E-06 

- 

3.58E-06 

- 

1.92E-07 

- 

1.85E-07 

- 


0.113 

8.56E-07 

3.37 

1.55E-06 

3.56 

6.58E-08 

4.56 

6.34E-08 

4.56 

2 

0.073 

2.27E-07 

3.06 

4.06E-07 

3.09 

5.65E-09 

5.65 

5.67E-09 

5.56 


0.038 

2.96E-08 

3.12 

5.30E-08 

3.12 

6.17E-10 

3.39 

5.97E-10 

3.45 


0.024 

6.87E-09 

3.15 

1.24E-08 

3.14 

7.78E-11 

4.47 

7.57E-11 

4.45 


0.143 

2.13E-08 

- 

3.00E-08 

- 

1.04E-08 

- 

9.98E-10 

- 


0.113 

7.16E-09 

4.64 

1.06E-08 

4.44 

3.33E-09 

4.86 

3.20E-10 

4.85 

3 

0.073 

1.32E-09 

3.89 

1.80E-09 

4.08 

1.89E-10 

6.61 

1.83E-11 

6.58 


0.038 

8.65E-11 

4.18 

1.20E-10 

4.14 

1.47E-11 

3.91 

1.40E-12 

3.95 


0.024 

1.25E-11 

4.17 

1.75E-11 

4.16 

3.52E-12 

3.09 

3.32E-13 

3.10 


Table 5 _ 

History of eonvergenee of the approximation in Example \3. 5 \ l) (smooth solution). 
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h 

\\eu ||i: 

error 

nt 

order 

lleqllb 

error 

nt 

order 

Ibslle 

error 

h 

order 

lleuHl: 

error 

int 

order 


0.143 

2.49E-03 

- 

2.20E-02 

- 

1.40E-03 

- 

2.53E-0 

- 


0.113 

1.81E-03 

1.35 

1.62E-02 

1.29 

7.08E-04 

2.92 

1.84E-03 

1.36 

0 

0.073 

l.llE-03 

1.11 

l.lOE-02 

0.90 

2.94E-04 

2.02 

1.12E-03 

1.14 


0.038 

5.75E-04 

1.01 

7.23E-03 

0.64 

1.63E-04 

0.91 

5.77E-04 

1.02 


0.024 

3.49E-04 

1.08 

5.73E-03 

0.50 

9.09E-05 

1.26 

3.50E-04 

1.08 


0.143 

4.04E-04 

- 

8.38E-03 

- 

4.29E-04 

- 

3.97E-04 

- 


0.113 

1.80E-04 

3.45 

5.60E-03 

1.72 

2.08E-04 

3.09 

1.89E-04 

3.15 

1 

0.073 

7.93E-05 

1.88 

3.38E-03 

1.16 

8.83E-05 

1.97 

8.07E-05 

1.96 


0.038 

4.52E-05 

0.86 

2.00E-03 

0.80 

4.82E-05 

0.93 

4.53E-05 

0.88 


0.024 

3.03E-05 

0.86 

1.63E-03 

0.45 

3.23E-05 

0.87 

3.03E-05 

0.87 


0.143 

1.55E-04 

- 

4.37E-03 

- 

1.77E-04 

- 

1.57E-04 

- 


0.113 

8.10E-05 

2.78 

3.02E-03 

1.58 

9.16E-05 

2.81 

8.12E-05 

2.82 

2 

0.073 

4.91E-05 

1.15 

1.72E-03 

1.30 

5.35E-05 

1.24 

4.91E-05 

1.16 


0.038 

2.70E-05 

0.92 

9.71E-04 

0.87 

2.87E-05 

0.95 

2.70E-05 

0.92 


0.024 

1.70E-05 

0.99 

8.32E-04 

0.33 

1.81E-05 

0.99 

1.70E-05 

0.99 


0.143 

7.94E-05 

- 

2.73E-03 

- 

9.13E-05 

- 

8.02E-05 

- 


0.113 

4.89E-05 

2.06 

1.84E-03 

1.68 

5.45E-05 

2.19 

4.92E-05 

2.08 

3 

0.073 

3.59E-05 

0.71 

1.09E-03 

1.21 

3.90E-05 

0.77 

3.60E-05 

0.72 


0.038 

1.79E-05 

1.07 

6.34E-04 

0.83 

1.90E-05 

1.10 

1.79E-05 

1.07 


0.024 

1.07E-05 

1.10 

5.15E-04 

0.45 

1.14E-05 

1.10 

1.08E-05 

1.10 


Table 6 

History of convergence of the approximation in Example \3.5\ )) (Non smooth solution). 


4. Elliptic interface problem. Let us now consider and interface S that di¬ 
vides the domain Q in two disjoint subdomains and as Figure 10 show. Then, 
problem ( |3.1[ ) becomes 


(4.1a) 

(4.1b) 

(4.1c) 

(4.1d) 

(4.1e) 

(4.1f) 


g|si 


-V ■ q = / in 11, 
q + K.\/u = 0 in 
u = go on 
q-n = qn on Fat, 
= sd on S, 

= sn on H. 


Here and are defined by 


:= {x — erP : x G H and e ^ 0}, 
:= {x — enP : x G H and e ^ 0}, 


where {j G {1,2}) is the unit outward normal unit vector of the subdomain , 
Sd G and sat G are prescribed jumps at the interface. At H we 

adopt the convention n := n^. 

For the sake of simplicity we assume dCl to be polygonal (if not, we apply the 
technique explained in previous section). However, the interface S is not necessarily 
piecewise flat. The numerical results provided in section ( |3.4| ) for a boundary value 
problem, suggested that the distance between the computational domain and the 
boundary should be of order O(h^) with a family of paths normal to the computational 
boundary. That is why we interpolate the interface S by piecewise linear segments. 
The computational interface, denoted by divides the computational domain 
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in two disjoint unions of elements and D^. {j G {1,2}) is defined as := 
{x — : X G T^h ^ ^ 0}? where is the unit outward normal vector of the 

computational domain D-^. 





Fig. 9. Approximation of the x-component of q Example \3.5\ (non-smooth solution). Columns: 
meshsize h = 0.143 and 0.024. Rows: Polynomial of degree k = 0,1 and 2. 


The main idea is to impose the jump of the scalar variable, denoted by 5^, on 
the computational interface On the other hand, the jump sn will be imposed at 


H by using the idea explained in Section 3.3 


Following the approach by m , the method HDG applied to the interface problem 
seeks an approximation {q^^Uh^Xh) G x Wh x such that 

(4.2a) ■y)Dh - {uh, ■ v)o^ + (A/i, v ■ n)dDh = 0, 

(4-2b) {w,Vh-qh)DH + {{Qh “ Qh)-n,w)dDH = 

(4-2c) • n, Ai)aD,,\(ruSh) = 0, 

(4-2d) {Xh,lJ.)rD = {9 d,Ij)td^ 

(4.2e) ■ n,//)r„ = (qjv, , 
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Fig. 10. Example domain Q divided in two regions and Q? by an interfaee E 


for all G V/i X Wh x M^. We still need to specify the jump of the normal 

component of gr at S. 


Here \h is a single-valued function, however the approximation of u must be 
double-valued on Then, similarly to m, we let Xh be the approximation of 

and consider as an approximation of Thus, we define 

(4.2f) Uh •= A/i + 


where defined on satisfies 


(4.3) 




1 on n ifdKn^hy^^ and K G D^, 
0 otherwise. 


To complete the method we define the numerical fiux as usual 


Qh'=Qh-^{^h-Uh)n on dDh. 

4.1. Approximation In order to define an approximation of 5^, we use the 
same transferring technique used for the Dirichlet data on a curved boundary (3.3). 
Let e G T^h such that e = dKiGdK 2 and, without loss of generality, assume that e lies 
completely inside of . We denote by the approximation (q^^Uh) restricted 

to the domain D-^. Now, for each cc G e, we observe that cr(cc) C Ai fl and then, 
according to the approximation given in (3.5f), 


(4.4) 


ulix)^ul{x)+ f K ^E^^{ql)-m, 

J a(x) 


where E^‘^{q\) is the standard extrapolation of q‘f^ to the whole 
(3.5d). Similarly, 


space defined in 


(4.5) 


ui{x)^ui{x)+ f K ^E^^qi)-m, 

J cr(x) 


In this case E^^{q\) = q\. 

Combining both equations, 

ul{x) - ul{x) « ul{x) - ul{x) + [ m- f K-'^E^^{ql) • m. 

J (j{x) J (j{x) 
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This expression suggest the following approximation 

(4.6) s’ltix) := Sd{x) + f ■ m - f K~^E^^ql) ■ m. 

J a{x) J a{x) 

4.2. Imposition of sn> For approximating sat we use the same idea that we 
applied for a Neumann boundary edge. For each interface edge e G we consider 
He C the part of H associated to e. We denote by Kl and K‘^ the element of 
and where e belongs. Then, we impose the following condition at the interface S: 

(4.7) (£;^=(g^)+£;^"(g^) = (sjv,//)e^ e M^CEe), 

where is defined similarly as in ( |3.7[ ). 

4.3. Numerical results: Interface problem. Finally, in this section we con¬ 
sider three numerical examples showing the performance of our technique in elliptic 
interface problems. Since the computational domains and do not exactly fit 
and we exclude from the computation of the errors the triangles intersecting the 
interface. Let the set of triangles whose faces are not interface edges. We measure 
the errors using the following norms || • 11 ^ 2 ( 6 ^) 


( 

l|e*llL2(£^) : = I hK\\PdU-Uh\\h^oK)\ • 

Vi^eD;,:i^nEh=0 / 

Example 4.1 (Elliptical-shaped domain). We first solve a Poisson equation in a 
the domain ft = (—1,1)^ divided by the elliptical interface E described by (x/0.8)^ + 
(^/0.4)^ = 1. We take K = I and 


u = 


cos(^) 

sin(7rx) sin(7r^) 


in ft^ 
in ft^ 


as exact solution. The source term, transmission and Dirichlet boundary conditions 
are obtained from this exact solution. 


In Tablej^the history of convergence for this example is displayed. Similarly to the 
examples involving Neumann boundary data, the order of convergence for u and q are 
optimal whereas the convergence of the numerical trace is suboptimal, i.e., 

Moreover, even though superconvergence of the postprocessed solution is lost, it 
provides a more accurate approximation of u. Figure pT] shows the approximation 
obtained with meshsizes of h = 0.072 and 0.018; and polynomial degree /c = 0, 1 and 
2 . 

Example 4.2 (Kidney-shaped domain). We now consider the same exact solu¬ 
tion as in previous example, but considering a kidney-shaped described by (2[(x + 
0.5)^ —X —0.5)^ — [(x + 0.5)^ +0.1 = 0. In despite of the changes of convexity 

of this geometry. Table shows similar accuracy on the approximations as the ones 
obtained in Example |4.1[ Figure shows the quality of the approximations of the 
scalar variable and its postprocessing obtained with a meshsize of h = 0.069 
and polynomial degree /c = 0, 1 and 2. As expected, provides a more accurate 
approximation of Uh without significantly increase the computational cost. 
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Ll2(D^) 


lledL 

L2(g 

h) 

lie.* II 


k 

h 

error 

order 

error 

order 

error 


order 

error 


order 


0.072 

2.37E - 01 

- 

3.53£; - 01 

- 

3.66£; - 

02 

- 

4.16£; - 

02 

- 


0.035 

1.22£; - 01 

0.94 

1.92£; - 01 

0.87 

2.01£; - 

02 

0.85 

2.16£; - 

02 

0.93 

0 

0.018 

5.97£; - 02 

1.03 

9.22E - 02 

1.04 

1.05£; - 

02 

0.93 

1.09£; - 

02 

0.98 


0.009 

2.98E - 02 

1.01 

4..66E - 02 

0.99 

5.60£; - 

03 

0.92 

5.67£; - 

03 

0.95 


0.004 

1.50E - 02 

1.00 

2.34.E - 02 

1.00 

2.84£; - 

03 

0.98 

2.86E - 

03 

0.99 


0.072 

1.98E - 02 

- 

4..20E - 02 

- 

1.75£; - 

03 

- 

2.2AE - 

03 

- 


0.035 

4..95E - 03 

1.97 

1.01£; - 02 

2.03 

1.63£; - 

04 

3.37 

2.72E - 

04 

3.00 

1 

0.018 

1.24.E - 03 

1.97 

2.36£; - 03 

2.07 

2.05£; - 

05 

2.96 

3.9QE - 

05 

3.12 


0.009 

3.12£; - 04 

2.01 

5.78E - 04 

2.05 

7.79E - 

06 

1.41 

8.92E - 

06 

1.95 


0.004 

7.85£; - 05 

2.00 

1.43£; - 04 

2.02 

1.24.E - 

06 

2.67 

1.22E - 

06 

2.73 


0.072 

1.44£; - 03 

- 

3.93£; - 03 

- 

1.25E - 

04 

- 

1.38E - 

04 

- 


0.035 

2.00£; - 04 

2.80 

5.24£; - 04 

2.86 

2.69E - 

05 

2.18 

2.78E - 

05 

2.47 

2 

0.018 

2.43£; - 05 

3.01 

5.99£; - 05 

3.10 

1.89E - 

06 

3.80 

1.9QE - 

06 

3.79 


0.009 

3.07E - 06 

3.01 

7.54£; - 06 

3.01 

3.00E - 

07 

2.67 

3mE - 

07 

2.72 


0.004 

3.92E - 07 

2.98 

9.54£; - 07 

2.99 

3.99E - 

08 

2.92 

A.OIE - 

08 

2.93 


0.072 

l.lOE - 04 

- 

3.02£; - 04 

- 

7mE - 

06 

- 

8mE - 

06 

- 


0.035 

7.78E - 06 

3.76 

2.16£; - 05 

3.75 

1.79E - 

07 

5.21 

3.9QE - 

07 

4.65 

3 

0.018 

4.A9E - 07 

4.08 

1.17£; - 06 

4.16 

6.33E - 

09 

4.78 

8.78E - 

09 

5.08 


0.009 

2.82E - 08 

4.02 

7.17£; - 08 

4.06 

4..60E - 

10 

3.81 

AmE - 

10 

4.18 


0.004 

1.80£; - 09 

3.98 

4.49£; - 09 

4.01 

1.93E - 

11 

4.59 

2.02E - 

11 

4.63 


Table 7 

History of convergence of the approximation in Example \4.1\ (elliptical-shaped) 




II®«IIl2(D^) 

l|eq \\l‘^(Dj^) 

lledh 

l2(£ 

h) 

lien* II 


k 

h 

error 

order 

error 

order 

error 


order 

error 


order 


0.069 

2.37E - 01 

- 

3.76E - 01 

- 

3.80£; - 

02 

- 

4.39£; - 

02 

- 


0.035 

1.23£; - 01 

0.97 

2.05£; - 01 

0.90 

1.97£; - 

02 

0.98 

2.12£; - 

02 

1.08 

0 

0.018 

6.07£; - 02 

1.03 

9.73E - 02 

1.09 

1.09£; - 

02 

0.86 

1.13£; - 

02 

0.92 


0.009 

3.01£; - 02 

1.01 

4.79£; - 02 

1.02 

5.69£; - 

03 

0.94 

b.77E - 

03 

0.97 


0.004 

1.51£; - 02 

1.00 

2.41£; - 02 

1.00 

2.89E - 

03 

0.99 

2.91E - 

03 

1.00 


0.069 

2.13£; - 02 

- 

4.35£; - 02 

- 

1.89E - 

03 

- 

2.b9E - 

03 

- 


0.035 

5.30£; - 03 

2.06 

1.08£; - 02 

2.05 

A.31E - 

04 

2.19 

A.9AE - 

04 

2.45 

1 

0.018 

1.33£; - 03 

2.03 

2.63£; - 03 

2.07 

imE - 

04 

2.10 

1.08£; - 

04 

2.23 


0.009 

3.28E - 04 

2.01 

6.24£; - 04 

2.07 

2.22E - 

05 

2.20 

2.27E - 

05 

2.24 


0.004 

8.28E - 05 

2.00 

1.56£; - 04 

2.01 

5.7AE - 

06 

1.97 

5.78E - 

06 

1.99 


0.069 

1.56£; - 03 

- 

3.72E - 03 

- 

1.39E - 

04 

- 

1.7AE - 

04 

- 


0.035 

2.02£; - 04 

3.02 

5.62£; - 04 

2.79 

1.78E - 

05 

3.04 

1.93E - 

05 

3.25 

2 

0.018 

2.58£; - 05 

3.02 

6.53£; - 05 

3.15 

2.53E - 

06 

2.86 

2.60E - 

06 

2.93 


0.009 

3.19£; - 06 

3.01 

7.76E - 06 

3.07 

A.17E - 

07 

2.60 

AA9E - 

07 

2.63 


0.004 

4.04£; - 07 

3.00 

9.80£; - 07 

3.01 

5.03E - 

08 

3.07 

5.0AE - 

08 

3.08 


0.069 

1.31£; - 04 

- 

3.50£; - 04 

- 

1.27E - 

05 

- 

lAOE - 

05 

- 


0.035 

7.96£; - 06 

4.13 

2.11£; - 05 

4.15 

9A2E - 

07 

3.84 

9.68E - 

07 

3.94 

3 

0.018 

4.92£; - 07 

4.08 

1.27£; - 06 

4.11 

A.93E - 

08 

4.61 

AAOE - 

08 

4.63 


0.009 

2.92E - 08 

4.06 

7.37E - 08 

4.10 

2.22E - 

09 

4.18 

2.23E - 

09 

4.19 


0.004 

1.87£; - 09 

3.99 

4.71£; - 09 

4.00 

lAAE - 

10 

3.98 

lAAE - 

10 

3.98 


Table 8 

History of convergence of the approximation in Example \4-2\ (kidney-shaped) 
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Fig. 11. Approximation of the scalar variable in Example \4. 1\ Columns: meshsize of h = 0.072 
and 0.018. Rows: Polynomial of degree k = 0^ 1 and 2. 


Example 4.3 (Thermal conductivity). Finally, considering the example provided 
by im, we simulate the heat distribution u at steady state, due to the heat source 
/, over the domain ft = ( — 1,1)^ divided by a circular interface of radius R = 0.5 
centered at the origin. The source term / and the thermal conductivity tensor are 
given by 

/(x, y) = —10{x^ + — I5x‘^{x‘^ + ~ ^5y^{x^ + 

and K. = Kjl in {j = 1, 2). The exact solution of this problem is 


u = 





in G 

■^) r ^ in ’ 

1^2 J 


and we consider ki = K 2 = 100. Dirichlet boundary condition on T is derived 
from the previous equation. In this case the jumps sd and sn are both equal to 
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zero. Table validates the optimal convergence rates of order for the heat 

distribution u and the flux q. Figure shows the approximated heat distribution 
considering meshes of size h = 0.072 and 0.018, and polynomials of degree /c = 0, 1 
and 2. 


L ^ 

i ^ 

u 

c 

o 

J 

c 

J 

c 


Fig. 12. Approximations Uh (left) and of the scalar variable u of Example J^..2 
meshsize h = 0.069. Rows: Polynomial of degree k = 0, 1 and 2. 


Columns: 


Remark 4.1. If the mesh is fine enough, the errors eu, eg and eu* can be 
computed in the entire computational domain since the quadrature points of a 
triangle K C will eventually lie in VtK This happens in all previous examples. In 
fact, we computed the errors \\eu\\L‘^(Dh): lkq||L2(D^) ctnd \\eu^\\L 2 (^Q^y Their behavior 
and magnitude are similar to ones displayed in the convergence tables. 

5. Conclusions. We have proposed a technique for high order approximation 
of boundary value problems in curved domains with mixed boundary conditions. We 
have provided numerical evidence suggesting that the technique performs properly if 
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the family of paths is normal to the computational boundary. A practical way to 
always satisfy this restriction is to define by interpolating F using only piecewise 
linear segments. Moreover, we have extend this technique to elliptic interface problems 
where the interface is not necessarily polygonal. We have presented numerical results 
indicating that the order of convergence of are optimal for the error of u and q if the 
interface is interpolated by piecewise linear segments. 



Fig. 13. Approximation of the scalar variable in Example \4-3\ (thermal conductivity). Columns: 
meshsize h = 0.072 and 0.018. Rows: Polynomial of degree k = 0^ 1 and 2. 
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Ll2(D^) 


lledL 

L2(g 

h) 

lie.* II 


k 

h 

error 

order 

error 

order 

error 


order 

error 


order 


0.072 

9.10E - 03 

- 

5.66E - 01 

- 

1.42E - 

03 

- 

1.63E - 

03 

- 


0.035 

7.30E - 03 

0.31 

3.09E - 01 

0.84 

7.58E - 

04 

0.88 

8.22E - 

04 

0.96 

0 

0.018 

4.56E - 03 

0.68 

1.46E - 01 

1.08 

4.04E - 

04 

0.91 

4.21E - 

04 

0.97 


0.009 

2.50E - 03 

0.87 

7.32E - 02 

1.01 

2.09E - 

04 

0.96 

2.13E - 

04 

0.99 


0.004 

1.28E - 03 

0.96 

3.60E - 02 

1.02 

1.06E - 

04 

0.97 

1.07E - 

04 

0.99 


0.072 

1.39E - 03 

- 

5.99E - 02 

- 

5.95E - 

05 

- 

1.55E - 

04 

- 


0.035 

4.51E - 04 

1.57 

1.45E - 02 

1.98 

1.23E - 

05 

2.20 

1.87E - 

05 

2.95 

1 

0.018 

1.36E - 04 

1.74 

3.36E - 03 

2.12 

2.34E - 

06 

2.39 

3.26E - 

06 

2.52 


0.009 

3.73E - 05 

1.88 

8.24E - 04 

2.04 

5.18E - 

07 

2.19 

6.20E - 

07 

2.41 


0.004 

9.42E - 06 

1.98 

2.06E - 04 

1.99 

6.66E - 

08 

2.95 

8.03E - 

08 

2.94 


0.072 

1.69E - 04 

- 

3.94E - 03 

- 

1.28E - 

05 

- 

2.07E - 

05 

- 


0.035 

2.33E - 05 

2.77 

4.62E - 04 

3.00 

9.34E - 

07 

3.66 

1.19E - 

06 

3.99 

2 

0.018 

3.40E - 06 

2.78 

5.36E - 05 

3.11 

1.05E - 

07 

3.16 

1.19E - 

07 

3.32 


0.009 

4.73E - 07 

2.87 

6.64E - 06 

3.04 

1.38E - 

08 

2.95 

1.45E - 

08 

3.07 


0.004 

5.94E - 08 

2.98 

7.84E - 07 

3.07 

1.39E - 

09 

3.30 

1.43E - 

09 

3.33 


0.072 

1.35E - 05 

- 

1.58E - 04 

- 

8.50E - 

07 

- 

1.24E - 

06 

- 


0.035 

8.36E - 07 

3.89 

7.36E - 06 

4.28 

1.80E - 

08 

5.39 

2.51E - 

08 

5.45 

3 

0.018 

5.64E - 08 

3.90 

4.05E - 07 

4.19 

1.03E - 

09 

4.13 

1.23E - 

09 

4.37 


0.009 

3.94E - 09 

3.87 

2.43E - 08 

4.09 

6.98E - 

11 

3.92 

7.39E - 

11 

4.08 


0.004 

2.49E - 10 

3.97 

1.41E - 09 

4.10 

2.02E - 

12 

5.10 

2.17E - 

12 

5.08 


Table 9 

History of convergence of the approximation in Example \4.3\ (thermal eonduetivity) 
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